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ABSTRACT 

We present ID, 2D, and 3D hydrodynamical simulations of core-collapse supernovae including a parameterized 
neutrino heating and cooling scheme in order to investigate the critical core neutrino luminosity (L cr it) required 
for explosion. In contrast to some previous works, we find that 3D simulations explode later than 2D simulations, 
and that L cr ; t at fixed mass accretion rate is somewhat higher in 3D than in 2D. We find, however, that in 2D 
L cr i t increases as the numerical resolution of the simulation increases. In contrast to some previous works, 
we argue that the average entropy of the gain region is in fact not a good indicator of explosion but is rather a 
reflection of the greater mass in the gain region in 2D. We compare our simulations to semi-analytic explosion 
criteria and examine the nature of the convective motions in 2D and 3D. We discuss the balance between 
neutrino-driven-buoyancy and drag forces. In particular, we show that the drag force will be proportional to 
a buoyant plume's surface area while the buoyant force is proportional to a plume's volume and, therefore, 
plumes with greater volume-to-surface area ratios will rise more quickly. We show that buoyant plumes in 2D 
are inherently larger, with greater volume-to-surface area ratios, than plumes in 3D. In the scenario that the 
supernova shock expansion is dominated by neutrino-driven buoyancy, this balance between buoyancy and 
drag forces may explain why 3D simulations explode later than 2D simulations and why L C rit increases with 
resolution. Finally, we provide a comparison of our results with other calculations in the literature. 
Subject headings: supernovae: general - hydrodynamics - neutrinos - stars: interiors 



1. INTRODUCTION 

Multidimensional phenomena play a critical role in the core- 
collapse supernova (CCSN) mechanism. Instabilities such as 
proto-neutron star convection (Epstein 1979; Burrows & Fryx- 
ell 1993), neutrino-driven convection (Herant et al. 1994; Bur- 
rows et al. 1995; Janka & Mueller 1996; Murphy et al. 2012), 
and the standing accretion shock instability (SASI, Blondin 
et al. 2003; Blondin & Mezzacappa 2006) have abetted explo- 
sions in 2D (Marek & Janka 2009; Muller et al. 2012a; Mueller 
et al. 2012) for progenitors that refuse to explode in spherical 
symmetry (Rampp & Janka 2000; Liebendorfer et al. 2001; 
Thompson et al. 2003). The additional degrees of freedom 
afforded by multiple dimensions can also increase the dwell 
time of matter in the post-shock region where the accreting 
matter experiences a net gain of neutrino energy resulting in an 
increased efficiency of neutrino heating. Taken together, these 
multidimensional effects lower the critical neutrino luminosity 
threshold for which explosions are obtained when comparing 
2D to ID (Murphy & Burrows 2008). In spherical symmetry, 
the radial stability of the supernova shock has been studied in 
detail (Burrows & Goshy 1993; Yamasaki & Yamada 2007; 
Pejcha & Thompson 2012; Fernandez 2012). The current state 
of the research into the CCSN mechanism has been reviewed 
by Janka (2012) and Burrows (2012). 

We lack a first principles understanding of how the afore- 
mentioned multidimensional effects result in a lowered critical 
threshold for explosion in 2D, however, it is attractive to extrap- 
olate the trend to 3D and postulate that the threshold for explo- 
sion will be reduced even further, perhaps permitting energetic 
explosions for the vast majority of progenitors. Nordhaus et al. 
(2010) found precisely this using a simplified prescription for 
neutrino heating and cooling and deleptonization. Nordhaus 
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et al. found that the critical luminosity in 3D was reduced 
by 15% - 25% as compared to 2D simulations for a 15 Mq 
progenitor. Using a similar parametric approach, Hanke et al. 
(2012) attempted to reproduce the results of Nordhaus et al. 
for both the 15 Mq progenitor and a 1 1.2 M Q progenitor. Con- 
trary to the findings of Nordhaus et al., Hanke et al. find that 
there is little difference between the critical luminosities in 2D 
and 3D, while recovering the result that the critical luminosity 
in 2D is significantly lower than that in ID (Murphy & Bur- 
rows 2008). While differences exist between the approaches 
of Nordhaus et al. and Hanke et al., the exact cause of their 
disparate results is unclear. Burrows et al. (2012) report that 
the results of Nordhaus et al. were beset by inaccuracies in the 
gravity solver in CASTRO that have since been corrected, and 
very recently Dolence et al. (2012) present new 3D CASTRO 
simulations that corroborate the basic result of Nordhaus et 
al.: 3D simulations explode earlier than 2D. The important 
question of whether the threshold for explosion is lower in 
three dimensions is, thus, still an open one. 

We are on the precipice of achieving 3D simulations of core- 
collapse supernovae with full spectral neutrino transport and 
adequate resolution, though a number of 3D CCSN simulations 
employing various approximations have already been accom- 
plished. The 3D simulations to-date have been facilitated by 
approximations, typically to the neutrino transport, or by low 
resolution. Here we briefly mention a non-exhaustive list of 
3D results relevant to CCSNe extant in the literature. The 3D 
calculations of Mueller & Janka (1997) and Khokhlov et al. 
(1999) did not include the effects of neutrinos and excised 
the PNS from the domain. A number of early 3D simulations 
utilized the smoothed-particle hydrodynamics approximation 
(Fryer & Warren 2002, 2004; Fryer & Young 2007), which has 
certain disadvantages over grid-based codes (see, e.g., Plewa 
2001; Agertz et al. 2007; McNally et al. 2012; Sijacki et al. 
2012). Studies of the SASI have made various approxima- 
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tions to achieve 3D simulations (Blondin & Mezzacappa 2007; 
Iwakami et al. 2008, 2009; Fernandez 2010). Simulations ne- 
glecting the effects of neutrinos and employing a simplified 
equation of state (EOS) have been used to study the ampli- 
fication of magnetic fields in 3D (Endeve et al. 2010, 2012). 
Hammer et al. (2010), using a neutrino "lightbulb" scheme 
(Scheck et al. 2006), followed the evolution of the supernova 
explosion all the way through the envelope of the progenitor 
in 3D and examined the asymmetric development of insta- 
bilities. Some studies have focussed on highly-magnetized 
progenitors with approximate (or no) treatments of neutrino 
effects (Kuroda & Umeda 2010; Winteler et al. 2012). The 
hydrodynamical kick imparted to the PNS has been studied in 
3D by Wongwathanarat et al. (2010, 2012). The rotational sta- 
bility of the PNS has been explored in 3D by Ott et al. (2005). 
Many 3D calculations have focussed on the emergent gravita- 
tional wave signal from CCSNe (Ott et al. 2007, 201 1, 2012b; 
Scheidegger et al. 2008, 2010b,a; Kotake et al. 2009; Miiller 
et al. 2012b). Some studies have focused also on the neutrino 
signal from 3D CCSNe (Lund et al. 2012; Ott et al. 2012b). 
Takiwaki et al. (2012) present low-resolution 3D simulations 
with spectral neutrino transport using the isotropic diffusion 
source approximation (IDSA, Liebendorfer et al. 2009). And 
recently fully-3D Boltzmann transport for neutrinos has been 
developed by Sumiyoshi & Yamada (2012). 

In this paper, we describe our multidimensional study of 
neutrino-driven CCSN explosions using a parameterization of 
the neutrino effects similar to that of Nordhaus et al. (2010) 
and Hanke et al. (2012). We find that the delay time until 
explosion for a given neutrino luminosity is greater in 3D than 
in 2D, i.e., L cr it is greater in 3D than in 2D. In Section 2 we 
describe our computational approach. In Section 3, we present 
our results. In Section 4 we discuss the dependence of the 
critical neutrino luminosity on dimensionality and resolution 
and suggest that our results can be understood by consider- 
ing the balance between buoyant and drag forces acting on 
neutrino-driven bubbles. We also examine the difference in the 
character of the shock motion between 2D and 3D. In Section 
5 we discuss the implications of our results and we conclude 
in Section 6. 

2. NUMERICAL METHOD 

Our numerical simulation approach is similar to that de- 
scribed by Couch (2012). We solve the Eulerian equations of 
hydrodynamics, 



dp 
dt 



V-O)=0, 



dpv 

Ik 
dpE 



(1) 

V • (pvv) + VP = -pV$, (2) 
V- [(pE + P)v] =pv- V$ + p(H-C), (3) 



where p is the mass density, v the velocity vector, P the pres- 
sure, $ the gravitational potential, E the total specific energy, 
H is the specific neutrino heating, and C is the specific neu- 
trino cooling. We use the directionally-unsplit hydrodynamics 
solver provided by the FLASH simulation framework (Dubey 
et al. 2009, , Lee et al., in prep.) to solve equations (1) - (3). 
We use third-order piecewise-parabolic spatial reconstruction 
(PPM, Colella & Woodward 1984) and a hybrid Riemann 
solver that uses the HLLE solver inside of shocks and the 
HLLC solver in smooth flow. We use a "hybrid" slope lim- 
iter that applies the monotonized central (mc) limiter to linear 



Figure 1. Volume rendering of entropy for values above 14 kg baryon -1 
for the 3D model with L„,52 = 1-7 at 850 ms post-bounce. The shock 
front is also volume-rendered in gray-scale. The blue sphere is a isodensity 
contour marking the edge of the proto-neutron star. The explosion develops 
in a non-symmetric manner with the PNS recoiling in the direction opposite 
the dominant direction of the explosion. The high-entropy buoyant plumes 
display a great deal of small-scale structure in 3D. 



wave families and the more diffusive minmod limiter to non- 
linear, self-steepening wave families. We use a monopole 
approximation to calculate the self-gravity of the flow. To 
facilitate comparison with Nordhaus et al. (2010) and Hanke 
et al. (2012), we use the Shen et al. (1998) equation of state 
(EOS), as implemented by O'Connor & Ott (2010). 

Our approach for neutrino heating and cooling is that de- 
scribed by Murphy & Burrows (2008), the same approach 
used by Nordhaus et al. (2010) and Hanke et al. (2012). The 
neutrino heating and cooling are given by, 



H = 1.544 x 10 



211 



10 52 ergs" 1 J \4MeV 



100 km 



erg 
g • s. 



and 

C = 1.399 xlO 20 



T V 



2 McV 



(Y p + Y n )e- T »- 



erg 



g- s 



(4) 



, (5) 



where L Ue is the electron neutrino luminosity (note it is as- 
sumed that L Pe — L Ue ), T„ e is the electron neutrino temper- 
ature, r is the spherical radius, (Y p + Y n ) is the sum of the 
neutron and proton number fractions, r Ve is the electron neu- 
trino optical depth, and T is the matter temperature. In all 
of our simulations, T„ e is set to 4 MeV. We approximate the 
neutrino optical depth by a piecewise fit based only on the 
density (see Couch 2012). This avoids the need to calculate 
radial integrals for the optical depth and is justified because the 
factor e~ T "" is included in equations (4) & (5) only as a cutoff 
to the neutrino source terms at high-densities. Differences in 
the implementation of this cutoff result in different normal- 
izations for the critical luminosity curves (J. Murphy, private 
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communication). Hanke et al. (2012) chose to adjust the neu- 
trino opacities used so that their ID critical curves matched 
those of Nordhaus et al. We have not. 

We follow the approach proposed by Liebendorfer (2005) 
for following the evolution of the electron fraction, Y e . In this 
approach, calibrated with ID Boltzmann transport simulations, 
Y e is dependent only on density. This is strictly only applicable 
during the pre-bounce collapse phase, however, we continue 
to use the density-dependent electron fraction approach post- 
bounce, as done by Nordhaus et al. (2010). We have found that 
neglecting any changes in Y e post-bounce results in substan- 
tially earlier explosions for a given neutrino luminosity. We do 
not include the entropy changes due to deleptonization given 
in Liebendorfer (2005). 

We use ID spherical, 2D cylindrical, and 3D Cartesian 
geometries with adaptive mesh refinement (AMR) as imple- 
mented in FLASH via PARAMESH (v.4-dev, MacNeice et al. 
2000). For this study we use a fiducial resolution at the max- 
imum refinement level of 0.7 km in each direction. We limit 
the maximum refinement level with radius such that a pseudo- 
logarithmic radial grid spacing is obtained. Our refinement 
limiter takes the form 

Axj > rjr, (6) 

where Ax^ is the grid spacing in the i -direction at refinement 
level £, r is the spherical radius, and r\ is a parameter that sets 
the effective angular spacing. If equation (6) is not satisfied by 
a given AMR block, further refinement of that block is prohib- 
ited. For our fiducial resolution we set the finest grid spacing 
to 0.7 km and r\ = 1.25%, resulting in an effective "angular" 
resolution of 0?54. In ID, the simulated domain spans km 
to 5000 km, in 2D the domain is km to 5000 km in cylindri- 
cal radius, R, and -5000 km to 5000 km in z, and in 3D the 
domain is -5000 km to 5000 km in each Cartesian dimension. 
At the outer spatial limits of the domain, we set boundary con- 
ditions that apply power-law profiles to density and velocity 
that approximate the stellar envelope outside the domain. Such 
boundary conditions are critically important to the results of 
the present study as simple "outflow" boundary conditions 
overestimate the mass accretion rate at late times, altering the 
explosion time for near-critical luminosities. This is because 
"outflow" boundary conditions enforce a zero-gradient condi- 
tion for the flow variables which mimics a flat density (etc.) 
profile outside the simulation domain artificially enhancing the 
mass flux into the domain from the boundary. 

We use the 15 M Q progenitor of Woosley & Weaver (1995) 
in all of our simulations. 

3. RESULTS: EXPLOSION TIMES 

We have run a series of ID, 2D, and 3D simulations in which 
we varied the driving neutrino luminosity. We start in the pre- 
collapsed progenitor phase and follow the evolution through 
collapse, bounce, shock stagnation and eventual revival. Figure 
1 shows a volume rendering of entropy and the shock surface 
in a 3D simulation at 850 ms post-bounce. In Table 1 we 
give the explosion delay times for our series of simulations 
and Figure 2 shows the average shock radii as a function of 
time post-bounce for a number of our simulations. Figure 3 
shows the critical luminosity curves as functions of both post- 
bounce explosion time and mass accretion rate at explosion. 
We consider a model to have exploded once the average shock 
radius exceeds 400 km and does not subsequently fall back 
below this value (as in Nordhaus et al. 2010; Hanke et al. 



Table 1 

Explosion times and accretion rates at time of explosion. 



0.5 km 0.7 km 0.7 km 

Li/ e ' ^exp -^exp ^exp -^exp ^exp -^exp 

(10 52 erg/s) (ms) (M Q /s) (ms) (M /s) (ms) (M /s) 



ID 



2.0 
2.1 
2.2 



2.3 


943 


0.153 


822 


0.170 






2.4 


538 


0.221 


554 


0.221 






2.5 


380 


0.262 


389 


0.262 






2.7 


216 


0.310 


212 


0.310 






2.9 


200 


0.314 


197 


0.317 










2D 




2D 




3D 


1.7 


713 


0.190 


388 


0.260 


821 


0.175 


1.8 


490 


0.233 


309 


0.274 






1.9 


313 


0.278 


291 


0.284 


403 


0.261 


2.0 






263 


0.294 






2.1 


247 


0.298 


222 


0.313 


238 


0.302 



a Electron-neutrino luminosity. 

Time after bounce of onset of explosion. A symbol indicates that the model 
does not explode during the simulated period of evolution. 
c Mass accretion rate at onset of explosion. 




Figure 2. Average shock radii as a function of time relative to bounce for 
three neutrino luminosities in 2D and 3D. Also shown for comparison is the 
shock radius from the ID simulation with L„,52 = 2.3. Universally, the 
shock expands more rapidly in 2D than in 3D. Increasing the resolution in 2D 
delays explosion, as shown by the cyan curve. 



2012), though other metrics, such as reaching a critical value 
of the ratio of advection time to heating time in the gain region 
(e.g., Fernandez 2012) or satisfying the 'ante-sonic' condition 
(Pejcha & Thompson 2012) may be used (for a comparison 
of the difference between these metrics, see Dolence et al. 
2012). We find that the critical luminosity curve is lowered in 
multidimensional simulations as compared with spherically- 
symmetric simulations, consistent with all previous similar 
studies (Murphy & Burrows 2008; Nordhaus et al. 2010; Hanke 
et al. 2012; Couch 2012). When comparing 2D to 3D, however, 
we find interesting and heretofore unprecedented behavior: at 
our fiducial resolution the 2D simulations consistently explode 
earlier than 3D simulations at the same neutrino luminosity. 
Figure 2 shows that for a given neutrino luminosity the average 
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Figure 3. Critical neutrino luminosity curves in both mass accretion rate (left) and post-bounce time (right) for our set of ID, 2D, and 3D simulations with the 15 
Mq progenitor. We reproduce the common result that the critical curves are lower in 2D than in ID, however, we find that for our fiducial resolution of 0.7 km the 
3D critical curves are higher than the 2D curves. Thus, our results indicate that obtaining explosion is more difficult in 3D. Critical curves for increased resolution 
ID and 2D simulations are also shown. Increasing resolution in ID results in almost no difference in the curves, whereas in 2D the higher-resolution simulations 
explode later and the curves are much closer to the 3D counterparts. 



shock radius expands more quickly in 2D than 3D. 

In 2D, the explosion time for a given luminosity is sensitive 
to the grid resolution used. Increasing the finest grid resolu- 
tion to 0.5 km and reducing the radial refinement limiter, ?y, 
to 0.94% results in a 2D criticality curve much more simi- 
lar to the 3D curve at the fiducial resolution of 0.7 km (and 
r] = 1.25%). Increasing the resolution in ID simulations re- 
sults in almost no change in the explosion times as a function 
of neutrino luminosity. This very importantly indicates that 
the cause of the resolution dependence is connected to an in- 
trinsically multidimensional process. At present, we lack the 
necessary computational resources to carry out a resolution 
study in 3D, though our results certainly indicate the necessity 
of such a study. Hanke et al. (2012) do carry out a resolu- 
tion study including 3D and also find that the explosion times 
are very dependent on grid spacing. For simulations with 400 
unequally-spaced radial zones, Hanke et al. find that increasing 
the angular resolution of their spherical grid results in earlier 
explosions in 2D but later explosions in 3D. Considering sim- 
ulations with at least 600 radial zones, the dependence of the 
explosion times on resolution is not always consistent amongst 
their results; increasing the angular resolution delays explosion 
in 2D for some of their models. It is important to note that 
their fiducial resolution (3?) is substantial coarser than ours 
(~ 0?54) and the finest resolution they use in 3D (1?5) is still 
coarser than our resolution. 

4. INTERPRETATION AND ANALYSIS 

4. 1 . Shock Expansion Driven by Buoyancy 

So then, what is the explanation of our results? That is, 
why is it that our 2D simulations explode earlier than our 
3D? and why does increasing the resolution in 2D result in 
later explosions? As we will argue in the following sections, 
our results indicate that neutrino-driven buoyant convection 
is the dominant instability that encourages shock expansion 
for this progenitor, particularly in 3D. Similar arguments have 
been made recently by Burrows et al. (2012) and Murphy 
et al. (2012). In this picture, accreting gas in the gain region 
absorbs neutrino energy eventually becoming buoyant and rises 




Figure 4. Constant entropy contours for a value of 14 kg baryon - 1 for 
models with L l/} $2 = 1-7 at the respective times of explosion (see Table 1). 
The left panel shows the 3D data and the right shows the 2D data revolved 
about the symmetry axis. As discussed in the text, 2D simulations show 
buoyant plumes that have much smaller surface area-to- volume ratios than for 
3D. 

toward the shock where this buoyant energy is used to push 
the shock further out. The speed at which a plume rises will be 
determined by the competition between the plume's buoyancy, 
determined by the amount of neutrino energy it has absorbed, 
and the drag force from cold gas traveling downward through 
the gain region (see, e.g., Dolence et al. 2012). A plume's 
buoyancy will be proportional to its volume since the neutrino 
energy absorption rate will scale as the solid angle of the plume 
times its neutrino optical depth. The drag force pushing back 
against the buoyant plume will scale as the surface area of 
the plume. Therefore, smaller plumes with greater surface 
area-to-volume ratios will rise more slowly than larger plumes. 

To see this, consider a spherical buoyant bubble of radius 
Tb at a distance from the coordinate origin R^. The instanta- 
neous buoyant force on this bubble is provided by neutrino 
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Figure 5. Entropy pseudo color plots for 3D and 2D simulations with L„ 52 = 1.7 at 100 ms post-bounce (left) and time of explosion (right). For 3D, three 
orthogonal slice planes are shown. By 100 ms the shock structure in 3D is still very spherical and high-entropy buoyant plumes are just starting to reach the shock. 
In 2D at 100 ms, elongation along the symmetry axis is already evident, particularly in the southern hemisphere, and the convective structures are larger and 
more coherent. At explosion time, 388 ms and 821 ms in 2D and 3D, respectively, the character of the shock and the buoyant convection behind it is completely 
divergent between 2D and 3D. The 2D explosion is characteristically dipolar and dominated by a few, large buoyant plumes (or arguably 1=1 SASI). In 3D, the 
shock does not show a dominant low-order shape and the convective plumes show much more small scale structure. There are also a greater number of low-entropy 
down flows in 3D. These fundamental differences between 2D and 3D result in convective plumes that have much greater surface area-to- volume ratios in 3D. This 
results in a greater drag-to-buoyant force ratio which, for buoyancy-dominated shock expansion, leads to slower shock expansion in 3D. 



radiation and so must be equal to the neutrino radiation force: 
F v ~ a v pbr\L v j cR\, where <f v is an effective neutrino cross- 
section that includes any geometric constants, pi, is the density 
of the bubble, and c is the speed of light. The instantaneous 
drag force on the bubble is Fd ~ CdpbV 3 rf, where Cd is a 
drag coefficient that contains any geometric constants, and v is 
the bubble's velocity relative to the background accretion flow. 
The ratio of the buoyant force to the drag force on the bubble 
is then 

F„ ^ o v L v r h 
F d ~ C d v 2 cR 2 b ' 

Thus, this ratio increases with bubble size, r^, and larger bub- 
bles or plumes will rise faster than smaller. 

In 2D, the typical plume size is larger than in 3D for multiple 
reasons. First, the axial symmetry intrinsic to 2D cylindrical 
coordinates means that off-axis plumes are really rings. In 3D, 
such plume rings are unstable and will break up into many 
smaller plumes. In Figure 4 we demonstrate this difference in 
plume scale with constant entropy contours in 3D (left) and 2D 
revolved around the symmetry axis (right). In 2D the forced 
symmetry results in very large "3D" plumes whereas in 3D no 
such large scale plumes can exist. These contours are plotted 



at the respective explosion times in correspondence with the 
right panel of Figure 5. Second, the symmetry axis in 2D 
encourages the growth of large plumes along it, due either to 
the action of low-order modes of the SASI or that of buoyant 
convection. In Section 4.2 we show that the amplitudes of low- 
order spherical harmonic modes of the shock deformation are 
reduced in 3D as compared to the 2D case. Third, the "inverse" 
turbulent energy cascade in 2D (Kraichnan 1967) will pump 
energy to larger scales whereas in 3D the "forward" energy 
cascade will send energy to smaller scales (see Section 4.4). 

In Figure 5 we compare entropy slices between 3D and 2D 
simulations with L v 52 — 1-7. The left half of Figure 5 is 
at a time of 100 ms post-bounce and the right half is at the 
time of explosion, 821 ms for 3D and 388 ms for 2D. At 
100 ms, the 3D simulation shows a shock that is still nearly 
spherical and developing convection behind it. The largest 
convective plumes are just reaching the shock perturbing it's 
spherical structure and stochastically pushing it out in radius. 
The 2D simulation is similar, but the developing convection 
in the gain region is visibly more coherent and vortex-like. 
Again, in 2D these vortical convective cells truly represent 
three-dimensional rings. Also in 2D, the influence of the sym- 



6 



COUCH 



0.16 
I 0.12 
3 0.08 

> 0.04 

0.16 

I 0.12 

S 0.08 h 
A 

> 0.04 





0.08 




o 








0.06 






0.04 




ST 




: 


> 


0.02 






0.0 0.2 0.4 0.6 0.8 0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 

tpb H t pb [s] t pb [s] 

Figure 6. Normalized power in the first three spherical harmonic modes of the shock deformation in 2D (red) and 3D (blue) for three neutrino luminosities. The 
"power" is the squared spherical harmonic coefficient summed over all m's, so in 2D where the m's cancel P(i) = <x? , The amplitudes for 3D are generally 
smaller than for 2D. The data necessary to compute the shock spherical harmonics for the 3D model with L v ^2 =1.7 prior to t = 0.3 s was, mistakenly, not 
saved. 



metry axis is already apparent as the shock is becoming elon- 
gated along it, particularly along the southern axis. At the 
time of explosion, the 2D and 3D structures have diverged 
significantly. The 2D explosion shows the characteristic "bipo- 
lar" shock structure and is dominated by about three large 
buoyant plumes, one each along the axes and one, somewhat 
smaller buoyant plume north of the equator. The 3D structure 
is quite different: many smaller buoyant plumes exist and are 
more evenly distributed in solid angle and the shock, while cer- 
tainly not perfectly spherical, does not show large-amplitude 
low-order deformation as in the 2D case. The salient point 
is that Figures 4 and 5 clearly demonstrate that the surface 
area-to-volume ratio of the buoyant plumes is much higher in 
3D than in 2D. Thus, due to the greater amount of drag relative 
to buoyant force the plumes will rise more slowly in 3D than 
in 2D and will, therefore, encourage a slower growth of the 
average shock radius and later explosion times. 

4.2. Character of the Shock Motion 

Multidimensional phenomena such as neutrino-driven con- 
vection and the SASI naturally result in aspherical shock mo- 
tion in CCSNe. Linear analysis indicates that the most unstable 
modes of the SASI will be low-order, I sa 1 (Foghzzo et al. 
2007; Guilet & Foglizzo 2012). Three-dimensional simula- 
tions (Iwakami et al. 2008, 2009) and 2D R — <f> axisymmetric 
simulations (Blondin & Shaw 2007) tailored to study the SASI 
show that the SASI develops an to = 1 "spiral" mode that can 
be considered the superposition of two or more out-of -phase 
£ = 1 modes (Fernandez 2010). Neutrino-driven convection 
will excite much higher-^ modes of the shock motion and we 
have in the previous sections argued that such buoyant convec- 
tion is the dominant instability that encourages expansion of 
the shock in the present simulations. In this section we justify 
this assertion by analyzing the shock motion via spherical har- 



monic decomposition. Our approach is comparable to that of 
Burrows et al. (2012). 

We can decompose some scalar quantity, X, into spherical 
harmonic components with coefficients 



x(e,4>)Yz m (d,<t>)<m, 



(8) 



where the spherical harmonics are 



( V2NpRp (cos 9) cos m<t> to > 0, 
Y™ = J N$P$ (cos 9) to = 0, (9) 

sin \m\4> to < 0, 



and 



NT = 



121 + 1 {t-m)\ 
4tt (£ + m)V 



(10) 



The spherical harmonics are computed from the associated 
Legendre polynomials, Pp, and we use the physics convention 
for spherical coordinates: 9 is the polar angle and <f> is the 
azimuthal angle. 

In order to study the shock behavior, we decompose the 
shock surface R s (9, (f>) according to equation (8) and define a 
spherical harmonic "power": 



(id 



where we follow Burrows et al. (2012) in normalizing the 
coefficients by a factor ( - 1 ) I m I / y/4Tr(2£ + 1) such that a 0Q = 
(R s ), an = (x s ), ai_i = (y s ), and a w = (z s ). In 2D, all 
m 7^ components cancel so that P(£) = a\. We use an 
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Figure 7. Mass-averaged entropy in the gain region. The average entropies 
in 3D are systematically higher than in 2D, a result of the smaller mass in 
the gain region for 3D. We do not find that average gain region entropy is an 
indicator of proximity to explosion. Shown for comparison are the ID data 
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accurate shock finding algorithm in FLASH to track and flag 
zones that are within the shock (Balsara & Spicer 1999). 

We show in Figure 6 the first three spherical harmonic pow- 
ers of the shock surface, normalized by the average shock 
radius, for both 2D and 3D simulations at various neutrino 
luminosities as functions of time. We find that the amplitudes 
of the low-order harmonics are significantly reduced in 3D as 
compared to 2D, particularly prior to the onset of explosion 
(marked by vertical dotted lines in Figure 6). In 3D, the am- 
plitudes are especially small prior to accretion of the Fe/Si 
interface at around 125 ms. When this occurs, the higher en- 
tropy in the Si shell encourages rapid shock expansion (Fig. 
2) that excites substantial shock deformation. This transient 
high-amplitude spike fades away, however, after the entirety 
of the Fe/Si interface is accreted through the shock. We find 
a general trend that the later the explosion time, the higher 
the P(l, 2, 3) at explosion. For example, the 3D model with 
L„ ( 52 = 2.1 explodes with very small P(l, 2, 3) and only after 
explosion sets in do the amplitudes grow to larger values. 

Of note, and not evident in Figure 6, is that between 50 
and 100 ms we find a low-amplitude spiral mode of the shock 
deformation. This spiral mode is damped once neutrino-heated 
plumes reach and impinge upon the shock at around 100 ms 
(see Fig. 5). Also, we have computed the spherical harmonic 
powers of the shock deformation for the higher-resolution 2D 
simulations and find no significant difference from the fiducial 
resolution 2D simulations. 

4.3. Explosion Indicators 

A number of quantities have been identified as possible indi- 
cators or predictors of explosion in core-collapse supernovae. 
Nordhaus et al. (2010) argued that a higher average entropy 
in the gain region correlates with likelihood of explosion and 
found a clear dimensional hierarchy in this quantity: 3D > 2D 
> ID. Hanke et al. (2012), however, suggest that there exists 
only a small, possibly insignificant, dependence of the gain 
region average entropy on dimensionality. In Figure 7 we plot 
the mass-averaged entropy in the gain region for our 2D and 
3D simulations. We define the gain region as any part of the 
simulation domain that has a net positive neutrino heating as 
determined by the difference between equations 4 and 5. We 
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Figure 8. Mass in the gain region as a function of time. Gain region mass is 
higher in 2D than in 3D, resulting in a higher net heating rate (Fig. 9) and a 
lower mass-averaged entropy (Fig. 7). 




Figure 9. Net neutrino heating rate in the gain region. Because of the greater 
gain region mass in 2D, the net heating is also higher. The net heating rises 
dramatically at late times (thick blue line) because the average temperature in 
the gain region drops, reducing the neutrino cooling rate (eq. 5). 

find that there is a clear hierarchy between 2D and 3D with 
3D having significantly higher average entropies than 2D. The 
ID average entropy, however, does not fit the dimensional hier- 
archy. Our results point out that average entropy in the gain 
region is not a good indicator of explosion likelihood as we 
find that 2D explodes before 3D yet still has smaller average 
entropies. The average entropy in the gain region is rather a 
measure of the time-integrated specific energy absorbed. Thus, 
the lower average entropies in 2D reflect that the gain region 
in 2D contains a greater amount of mass than in 3D, as shown 
in Figure 8. This higher gain region mass is also reflected in 
a higher integrated heating rate (Fig. 9). A higher net heat- 
ing rate should, intuitively, result in a greater rate of shock 
expansion. 

Pejcha & Thompson (2012) suggest that instability to run- 
away shock expansion sets in when the maximum value of the 
sound speed squared to the escape velocity squared reaches 
a critical value, the so-called 'antesonic' condition. For 
an isothermal equation of state, the antesonic condition is 
max(Cs/i)p SC ) = 0.19, where c s is the adiabatic sound speed 
and Vosc is the escape speed. For a polytropic equation of 
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Figure 10. Maximum value of the ratio of the squared adiabatic sound speed 
to the squared escape velocity times adiabatic index, the 'antesonic' condition. 
We find that the antesonic condition is a good predictor of proximity to 
explosion, although our data indicate that the critical value of ~ 0.2 suggested 
by Pejcha & Thompson (2012) is somewhat low. We also plot the data for the 
high-resolution 2D simulations. The antesonic condition is a better predictor 
of the delayed explosion time for the high-resolution simulations than any of 
the other integral quantities plotted in Figs. 7-9. We also show the critical 
ratio for the ID model with L„,52 = 2.3 (black dotted line). The ID data 
look very similar to the 3D data for -L„,52 = 1.7 and peak at explosion time. 

state, the antesonic condition becomes max(c^/fc SC ) = 0.19I\ 
where T is the polytropic index. For the completely general 
EOS we use, we modify the polytropic antesonic condition to 
be 

max( C 2/ Uc 2 sc7c ) = 0.19, (12) 

where 7 C is now the varying adiabatic index given by the EOS. 
In Figure 10 we plot the left hand side of equation (12) as a 
function of time for our simulations. We find that the antesonic 
condition is a good indicator of the beginning of accelerated 
shock expansion, although a critical value of 0.2 may be a 
little low (see also Hanke et al. 2012; Dolence et al. 2012). We 
see that an antesonic value of about 0.3 corresponds very well 
with our definition of the explosion time, when the average 
shock radius exceeds 400 km. Also, the antesonic value for 
our 2D simulations exceeds the critical value before the 3D 
simulations, echoing our result that 2D explodes before 3D. 

In Figure 10 we also plot the antesonic value for 2D simu- 
lations with higher resolution. Interestingly, of the quantities 
shown in Figures 7-10 the antesonic condition shows the 
most noticeable differences between the 0.7 km and 0.5 km 
resolution 2D simulations. The difference is most notable for 
the simulations with = 1.7 where we see that based on 
the antesonic condition we would expect the higher resolution 
simulation to explode later as, in fact, it does. 

4.4. Non-Radial Motion and Turbulence 

Non-radial motion, particularly on large scales, has been 
suggested as a primary factor resulting in easier explosions in 
multidimensional simulations as compared with spherically- 
symmetric calculations (Hanke et al. 2012). Such motion 
increases the matter dwell times in the gain region and, thus, 
the net neutrino heating rate. A measure of non-radial motion 
is the gain region kinetic energy in transverse, (9/</>direction, 
motion, which we plot in Figure 1 1 . The kinetic energy of 
transverse motion in the gain region is seen to corollate with 
increased neutrino luminosity. The transverse kinetic energy 
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Figure 11. Non-radial components of the kinetic energy in the gain region. 
The transverse kinetic energy is higher in 2D than in 3D, but this is again a 
reflection of the greater mass in the gain region for 2D. 




Figure 12. Energy spectra of the ^-direction kinetic energy, a proxy for the 
turbulent kinetic energy, in spherical harmonic basis calculated in a ~ 10 km 
wide shell centered on 150 km. We plot spectra for both 2D and 3D, and 
the spectra are averaged over 10 ms centered on the times indicated. For 3D 
we find an inertial range from about I = 10 — 70 over which the spectrum 
is roughly consistent with a £ -5 / 3 power-law. Dissipation sets in around 
I ~ 70 where the spectrum falls off more steeply than £ -5 / 3 . In all cases, 
the 2D spectra show more energy at i = 1 than the 3D spectrum, reflective 
of the inverse energy cascade in 2D. At around the driving scale of I ~ 10, 
the 2D spectra follow a t~ 3 power-law consistent with the forward enstrophy 
cascade. The dissipation scale for the higher resolution 2D simulations is at 
noticeably larger I. 

is also typically higher in 2D than in 3D. 

Turbulence has been suggested to play an important role 
in the supernova shock expansion in multidimensional simu- 
lations (Murphy & Meakin 201 1; Murphy et al. 2012). The 
nature of turbulence between 2D and 3D, however, is funda- 
mentally different. The best-known example of this funda- 
mental difference is the so-called 'inverse' energy cascade in 
2D: energy is transported to large scales in 2D whereas in 3D 
energy is transported to smaller scales (Kraichnan 1967). The 
characteristic power-law slope of the energy cascade is -5/3 
in either the spherical harmonic mode, I, or wavenumber, k. 
Also in 2D, enstrophy, a quantity proportional to the squared 
vorticity, is transported to smaller scales in a so-called forward 
cascade with a characteristic power-law index of -3. 
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In order to study the character of the turbulence in our simu- 
lations we have computed the ^-direction kinetic energy spec- 
trum in spherical harmonic basis where the energy in spherical 
harmonic mode £ is 

i 

E(i) = £ 4 m , (13) 

and the coefficients ag m are given by equation (8) with 
X(9, 4>) = i//o(r, 8, cf>)ve(r, 8, <fi). Transverse kinetic energy 
spectra for our multidimensional simulations are shown in 
Figure 12 where we restrict the evaluation to a ^10 km shell 
centered at r = 150 km. We also average the spectra over 
10 ms centered on the respective labeled times. The differing 
character of turbulence between 2D and 3D is evident. Begin- 
ning around £ — 10, which we identify as the driving scale 
of the turbulence, the 3D spectrum is consistent with a £~ 5 / 3 
power-law and the 2D spectra roughly follow a £~ 3 power-law, 
consistent with our expectations based on turbulence theory 
(Kolmogorov 1941; Kraichnan 1967). Also evident is that there 
is much more energy at small scales (large ts) in 3D than in 
2D, while 2D has more energy at large scales {1=1 — 3). 
In a real sense, the different natures of these energy spectra 
reflect the different characteristic buoyant plume sizes seen 
in our simulations: plumes are smaller and more numerous 
in 3D as compared with 2D. This is then directly related to 
our argument that the smaller plumes in 3D will experience 
a greater drag-to-buoyant force, slowing their average ascent 
relative to the 2D case, resulting in a less-rapid average shock 
expansion in 3D. 

Figure 12 also shows energy spectra for 2D at different times 
and at different resolutions. Note that for all the spectra plotted 
the average shock radius is approximately the same at the time 
shown (250-275 km, see Fig. 2). Considering the evolution of 
the 2D energy spectra with time, the inverse energy cascade is 
evident in the increase in energy by an order of magnitude at 
£ = 1 while the spectrum in the inertial range, £ > 10, remains 
relatively unchanged in the later time. The lower-resolution 
(0.7 km) 2D spectrum shows two important distinctions. First, 
dissipation sets in at smaller £, consistent with the assumption 
that grid spacing determines the dissipation scale. Second, 



at 200 ms there is substantially more energy on large scales 
(£ = 1,2) than for the 0.5 km resolution case. 

It is of curious, perhaps surprising, significance that the en- 
ergy spectra seem to match the expected power-law slopes 
so well. The power-law scaling in the inertial ranges of -5/3 
for 3D and -3 for 2D are based on the assumption of fully- 
developed isotropic turbulence (Kolmogorov 1941; Kraichnan 
1967). The turbulence in our simulations is not isotropic but 
convective. There is no model of buoyant convective turbu- 
lence that would suggest any other power-law scalings, but 
that they seem to be the same is an interesting coincidence. 

4.5. Resolution Dependence 

As discussed in the previous sections, we find that increasing 
the resolution in 2D simulations results in later explosion times. 
This is an expected result of equation (7) that smaller buoyant 
plumes rise more slowly than large buoyant plumes due to the 
greater amount of drag force relative to buoyant force. The 
different convective plume size is evident by the differences in 
the energy spectra in Figure 12, but also by visual inspection of 
the 2D data. In Figure 13 we show a side-by-side comparison 
of the 2D simulations with different resolution at two different 
post-bounce times. What is evident in this comparison is that, 
even at 200 ms, the higher-resolution simulation shows a larger 
number of low-entropy down flows and deep penetration of 
these down flows. The larger number of down flows breaks 
up the rising buoyant plumes while their deeper penetration 
further increases the surface are of the plumes and, therefore, 
the drag force felt by the plume. This is especially evident at 
explosion time (right half of Fig. 13). 

This is simple to understand by considering the primary 
instabilities that will contribute to break-up of the buoyant 
plumes. The plumes are inherently Rayleigh-Taylor unstable 2 
and will also develop parasitic Kelvin-Helmholtz instabilities 
as they rise through the post-shock accretion flow. In the in- 
viscid limit, the growth rate for both of these instabilities is 
inversely proportional to the grid scale (Chandrasekhar 1961; 
Youngs 1984), thus the higher-resolution simulation experi- 

2 In a sense, it is the Rayleigh-Taylor instability that drives the convection 
in the first place. 
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ences faster growth of these instabilities that impede the rising 
plumes, slowing overall shock expansion. 

Hanke et al. (2012) report that increasing the resolution in 
3D results in delayed explosions while increasing resolution 
in 2D results in earlier explosions. Their 2D resolution study 
results are contradictory to our simulation results. Hanke 
et al.'s conclusion that increasing resolution in 2D results in 
earlier explosions is based on their simulations using only 
400 radial zones. When considering their results with 600 or 
800 radial zones, however, the dependence on resolution is 
less clear; for some sets of 2D models with more than 400 
radial zones increasing the angular resolution results in later 
explosions, as we have found. Hanke et al. suggest that this 
behavior with increased radial resolution is due to a artificial 
density peak in the cooling region that grows with increasing 
radial resolution, enhancing total cooling. We do not see 
such an artificial density peak in our simulations. It is also 
worth noting that our fiducial "angular" resolution of 0?54 is 
substantially finer than the fiducial resolution considered by 
Hanke et al. (3?) and comparable to the finest resolution they 
use in 2D (0?5). Nordhaus et al. (2010) and Dolence et al. 
(2012) use a fiducial resolution greater than ours (0.5 km) but 
do not conduct resolution studies so it is unknown how their 
results would vary with resolution. 

5. DISCUSSION 

The results of previous studies similar to ours have been 
contradictory. The Princeton group has found that 3D simula- 
tions explode sooner than 2D (Nordhaus et al. 2010; Dolence 
et al. 2012) while the Garching group finds little difference in 
the explosion times between 3D and 2D (Hanke et al. 2012). 
The initial aspiration of our study was to cast a tie-breaking 
vote in this discussion but instead it has raised new questions 
by yielding a third result: our 3D simulations explode later 
than our 2D simulations. The source of the disparity in the 
results from the different groups is still unclear. There are sig- 
nificant differences in numerical approach between the Princ- 
ton and Garching groups. The Princeton group uses the new 
code CASTRO (Almgren et al. 2010) which implements a 
directionally-unsplit piecewise parabolic method with an HLL- 
type Riemann solver. CASTRO uses patch-based adaptive 
mesh refinement (c.f., Berger & Oliger 1984; Berger & Colella 
1989) as provided by the Box lib framework. The Princeton 
simulations are run in ID spherical, 2D cylindrical, and 3D 
Cartesian coordinates. The Garching group uses the venera- 
ble hydrodynamics code PROMETHEUS (Fryxell et al. 1991) 
which implements directionally-split PPM with a 'two shock' 
Riemann solver (Colella & Glaz 1985) in smooth flow and the 
HLLE Riemann solver in shocks. PROMETHEUS does not 
use AMR but instead relies on ID, 2D, and 3D spherical coor- 
dinates with non-equidistant radial spacing. Though perhaps 
the most important differences between the approaches of the 
Princeton and Garching groups are those dealing with approx- 
imations to the neutrino physics. Nordhaus et al. (2010) and 
Dolence et al. (2012) follow closely the approach of Murphy 
& Burrows (2008): deleptonization is approximated using the 
density-dependent parameterization of Liebendorfer (2005), 
both pre- and post-bounce, and post-bounce neutrino heating 
and cooling are considered locally based on rates derived by 
Janka (2001). Heating and cooling are shut off at high density 
by adding a e~ T " term to the rates. Hanke et al. (2012) follow 
a very similar approach except for the following. They follow 
collapse and bounce to 15 ms post-bounce with full neutrino 
transport in ID and do not employ the parameterization of 



Liebendorfer (2005) at all. The subsequent multidimensional 
evolution uses the same heating and cooling rates as in Nord- 
haus et al. (2010) but the neutrino optical depth is computed 
by integrating the appropriate opacity in radius whereas the 
Princeton group uses a density-dependent parameterization for 
the optical depth (J. Murphy, private communication, see also 
Couch 2012). In their multidimensional post-bounce calcula- 
tions, the Garching group follows the evolution of the inner 
core in ID to avoid issues with grid convergence. 

Our approach is very similar to those of Nordhaus et al. 
and Dolence et al.: directionally-unsplit PPM in ID spher- 
ical, 2D cylindrical, and 3D Cartesian with AMR. We use 
the Liebendorfer parameterized deleptonization pre- and post- 
bounce and approximate the neutrino optical depth with a 
density-dependent piecewise fit, although we scale optical 
depth differently than the Princeton group resulting in a differ- 
ent normalization of the critical curves (see also the discussion 
in Hanke et al.). Notable differences in our scheme and that 
of the Princeton group are we use a oct-tree, block-structured 
AMR package that does not sub-cycle in time, i.e. all refine- 
ment levels advance with the same time step size. CASTRO 
uses adaptive time refinement allowing coarser resolution lev- 
els to advance with larger time steps. And we use a "hybrid" 
Riemann solver that uses the HLLC solver in smooth flows 
and the HLLE solver in shocks. Another possibly significant 
difference amongst all three groups is in the EOS implementa- 
tion. While each group is using the Shen et al. high-density 
baryonic EOS for the referenced simulations, the construction 
of the tables actually used is no doubt somewhat different. 
There could also be important differences in the details of the 
implementations of the monopole gravity solver amongst the 
different codes. What is clearly mandated by the disparity in 
the results is a rigorous code-to-code comparison. 

Our results, particularly in 3D, support the conjecture that, 
for this progenitor and treatment of neutrino physics, the SASI 
is subdominant to neutrino-driven convection in advancing the 
average shock radius (Burrows et al. 2012; Murphy et al. 2012). 
We are, however, very hesitant to extrapolate this conclusion 
to all possible progenitors. For instance, Mueller et al. (2012) 
show clear evidence for a strong SASI in the 2D explosion of 
a 27 Mq progenitor using conformally-flat general relativistic 
dynamics and full neutrino transport. Although recently Ott 
et al. (2012a) have simulated the same 27 Mq progenitor with 
full 3D GR and a multispecies neutrino leakage scheme and 
found that neutrino-driven convection becomes the dominant 
instability in exploding models. 

In the scenario that neutrino-driven convection dominates, 
the shock expansion can be roughly described by buoyant 
plumes plumes rising against the drag force exerted by the 
post-shock accretion flow (Dolence et al. 2012). The buoyant 
force exerted by the neutrino radiation on a plume will be pro- 
portional to the plumes subtended solid angle times the optical 
depth of the plume multiplied by the neutrino luminosity, i.e., 
the plume's volume. The drag force exerted by the accretion 
down flows will be proportional to the plume's surface area. A 
natural result of this model is that a single plume with volume 
V will rise more quickly than two plumes with volumes V/2 
due to the greater amount of total surface area, and hence drag 
force, for the two plumes [see equation (7)]. As discussed at 
length in Section 4, 3D simulations develop smaller, more nu- 
merous buoyant plumes than 2D simulations. Thus the balance 
between buoyancy and drag explains why the 3D models ex- 
plode later than the 2D models. This picture also accounts for 
the 2D resolution dependence as higher-resolution will result 
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in plumes that have greater surface area-to-volume ratios. Sim- 
ilar behavior of delayed or precluded explosion with increased 
resolution has been seen in more realistic 2D simulations using 
full neutrino transport (Marek & Janka 2009). 

6. CONCLUSIONS 

We have conducted a ID, 2D, and 3D parameter study of 
neutrino-driven core-collapse supernovae designed to explore 
the difference that 3D makes on the explosion characteristics, 
particularly time of explosion relative to ID and 2D. We find, 
as have a number of previous studies, that the so-called 'crit- 
ical curve' in the neutrino luminosity-mass accretion rate at 
explosion time plane is significantly lowered in multiple di- 
mensions relative to the spherical-symmetric case. A novel 
result of our study is we find that 3D explosions occur later 
than 2D explosions at the same neutrino luminosity, i.e., the 
3D critical curve is higher than the 2D critical curve. We find 
that the 2D results are resolution-dependent: increasing the 
resolution in 2D delays explosion pushing the high-resolution 
2D critical curve very near to the 3D critical curve (Fig. 3). 

We suggest that our results can be explained by the compe- 
tition between buoyancy and drag. The shock expansion in 
our simulations is dominated by the action of neutrino-driven 
buoyant convection (see also Burrows et al. 2012; Murphy 
et al. 2012). In this case, the shock expansion can be fit by 
the motion of a buoyant plume rising through the shocked 
accretion flow (Dolence et al. 2012). The buoyant force acting 
on the plume, provided by the absorbed neutrino radiation, 
is proportional to a plume's volume whereas the drag force 
resulting from the downward-flowing accretion is proportional 
to a plume's area. Thus, a plume's ascension velocity increases 
with increasing volume-to-surface area ratio. Rapid shock ex- 
pansion, therefore, is best abetted by plumes that subtend large 
solid angles. We have shown that 3D simulations naturally 
result in many more, smaller-solid angle plumes than compa- 
rable 2D simulations. We posit that this, then, is why our 3D 
simulations explode later than our 2D simulations and why 
higher-resolution 2D simulations also result in later explosions. 

We examined several differences between our 2D and 3D 
simulations. In Section 4.2 we explored the character of the 
shock motion by calculating the first few spherical harmonic 
powers of the shock deformation in 2D and 3D. We find that the 
amplitudes of the low -order, I = 1,2, modes of the shock mo- 
tion that are often associated with the SASI are much reduced 
in 3D relative to 2D. This result is in qualitative agreement 
with Burrows et al. (2012). 

When considering possible indicators of explosion, we find 
that the "antesonic" condition of Pejcha & Thompson (2012) 
correlates well with proximity to explosion. We find that the 
mass-averaged entropy in the gain region is typically higher 
in 3D than in 2D, but our results indicate that is is not a good 
indicator of proximity to explosion, as it was suggested to be 
by Nordhaus et al. (2010). Instead this higher average entropy, 
which is proportional to the time-integrated neutrino heating 
per mass, only reflects that there is less mass in the gain region 
in 3D as compared to 2D. The greater gain region mass in 
2D also results in a greater net heating rate, which also aids 
explosion relative to 3D. 

The character of the non-radial motion and turbulence be- 
tween our 2D and 3D simulations is distinct. In 3D, the for- 
ward energy cascade transports energy to smaller and smaller 
scales until dissipation do to finite grid size sets-in. In 2D, 
the forward enstrophy cascade results in a much less efficient 
transport of energy to small scales while the inverse energy 



cascade actually transports energy from the driving scales to 
larger scales. This behavior is conducive to explosion as it 
encourages both the growth of the low-order SASI and the 
development of large buoyant plumes. 

The cause of the differences in results among the groups 
attempting similar multidimensional parameter studies (e.g., 
Nordhaus et al. 2010; Hanke et al. 2012; Dolence et al. 2012, 
this work) remain to be explained. The simplified physics 
employed by such studies would make a detailed code-to-code 
comparison effort far more straight-forward than a compar- 
ison that included neutrino transport. Such a code-to-code 
comparison will hopefully be seen as a high priority and be 
accomplished in the near future. 

The approximations we employ are admittedly crude. They 
have tremendous merit, however, in that they facilitate 3D 
parameter studies of CCSNe. While certain conclusions about 
CCSNe based on studies such as the present should be made 
with caution, the 3D results appearing in the literature to- 
date make one undeniable point: 3D core-collapse supernovae 
are fundamentally and dramatically different than 2D core- 
collapse supernovae. The absence of forced-symmetry in 3D 
makes an enormous impact on the character of the shock mo- 
tion and the development of neutrino-driven convection and 
turbulence. Our results indicate, however, that 3D alone may 
not be the key to successful, robust explosions in more realistic 
simulations. Still, the enormous difference between 2D and 
3D CCSN simulations emphasizes and underlines the need for 
fully 3D simulations by this community. 
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